Inhomogeneous relaxation in vibrated granular media: consolidation waves 
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We investigate the spatial dependence of the density of vi- 
brated granular beds, using simulations based on a hybrid 
Monte Carlo algorithm. We find that the initial consolida- 
tion is typically inhomogeneous, both in the presence of a 
constant shaking intensity and when the granular bed is sub- 
mitted to 'annealed shaking'. We also present a theoretical 
model which explains such inhomogeneous relaxation in terms 
of a 'consolidation wave', in good qualitative agreement with 
our simulations. Our results are also in qualitative agreement 
with recent experiments. 



I. INTRODUCTION 



In recent years the issue of compaction in granular media 
has attracted experimental ( [|lj, [^, [^]) as well as theoretical 



X LLj' LLj' fl)' interest in the physics community. In 

particular the important issue of compactivity [|l^, as an ef- 
fective temperature forpowders in the quasi- static regime, has 
intrigued researchers |3| who have sought, as a consequence, 
to investigate further the implications of such a definition. 

While this effective temperature was originally defined in 
terms of thermodynamic analogies for a static powder, some 
inconsistencies have motivated work on other possible 
definitions. A Langevin dynamics approach, used to model 
vibrated powders [O, was able to introduce the compactiv- 
ity as the natural interpolation of the conventional 'granular 
temperature' ||l2| as the intensity of vibration was decreased; 
this led to its interpretation as the effective temperature cor- 
responding to the slow dynamical modes of the system. This 
interpretation of the compactivity was recently reformulated 
at a more microscopic level, in terms of a p-spin model of 
granular dynamics | |l3[ . 

As the compactivity parameter is closely related to den- 
sity variations recent experimental investigations have 
focused on spatiotemporal density fluctuations in a shaken 
powder [0]. Recent molecular dynamics simulations have 
also highlighted this issue by focusing on the necessarily inho- 
mogeneous density fluctuations which arise when a granular 
system is forced by the motion of its boundaries. In this pa- 
per, we examine two aspects of spatial inhomogeneities: 

• First, we examine spatial density fluctuations in the 



presence of a constant shaking intensity, and find def- 
inite evidence of inhomogeneous relaxation at least at 
short times. Additionally, as in experiment, our 
simulations suggest that the degree of inhomogeneity is 
less at higher intensities of shaking. We put our obser- 
vations into perspective by suggesting that the density 
of a vibrated granular bed initially relaxes by the prop- 
agation of a 'consolidation wave'; this behaviour can 
be represented by a theoretical model that is in good 
qualitative agreement with our simulations. 

Second, we examine the response of the granular bed to 
an 'annealed shaking' procedure ||2|; the consolidation 
along the socalled reversible branch of the compaction 
curve shows evidence of strong spatial inhomogeneities, 
which we demonstrate via the configurational histories 
of the particles. Compaction proceeds from the walls, 
and propagates Uke a consolidation wave; this is rem- 
iniscent of earlier studies of transitions to crystallinity 
[l6| in shaken granular beds, and has also been ob- 
served in a recent experiment |l7| on the appearance of 
'nematic' ordering in an assembly of rods. 



II. PROCEDURE; THE HYBRID MONTE CARLO 
MODEL 

Our simulations use uniform hard spheres, subjected to 
non-sequential reorganizations which represent the effect of 
shaking. A variable shaking amplitude A is parametrised in 
units of the particle diameter; thus for example, A = 1.0 
means that shaken particles are able to move longitudinally 
and laterally by, 1 particle diameter (subject to volume ex- 
clusions) during a shake cycle. The details of the shaking 
algorithm have been discussed elsewhere ( jlj]). Briefly, 
one cycle of vibration of the granular assembly (corresponding 
to one timestep of our simulation) is modelled by: 

1. a vertical dilation of the granular bed, in proportion to 
the shaking amplitude A 

2. a stochastic rearrangement of the individual particles 
in an external field and with available free volume pro- 
portional to the shaking amplitude 
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3. and finally a cooperative recompression of the assem- 
bly as each grain lands on the substrate alone or with 
neighbours (in the latter case, arches could form). 

This algorithm has both static and dynamic disorder. The 
former is incorporated via the non-lattice-based nature of the 
simulations, while the latter is a consequence of the stochastic 
step, which allows particles the freedom to move in planes 
perpendicular to the direction of vibration, in proportion to 
the driving force. Collective structures, such as arches, are 
also allowed to form and propagate, because of the inclusion 
of the cooperative recompression step. In the absence of this 
step, particles would land sequentially, and the configurations 
would include only simple three-point stability conditions. 

III. RESULTS ON SPATIAL DEPENDENCE OF 
PACKING FRACTION 

In this section, we explore first numerically and then ana- 
lytically, the spatial dependence of the packing fraction of a 
shaken powder. We were guided in our investigations by the 
knowledge that typically, experiments rely on indirect knowl- 
edge of spatial variations of material properties, since the use 
of 3d techniques, such as MRI to make quantitative mea- 
surements, is still not very widespread. Consequently, we look 
at macroscopic quantities such as the average packing fraction 

< (j) > and the mean particle height Z(t), and extract infor- 
mation on spatial inhomogeneities of < >, indirectly . In 
the analytical section we discuss related issues, once again 
from a macroscopic point of view. 

A. Numerical simulations 

In Fig. |l| we have plotted the mean packing fraction 

< > against time t for two different amplitudes A, where 
the angular brackets indicate an ensemble average over dis- 
tinct realizations (independent sequences of pseudo-random 
numbers) of the noise. 



FIG. 1. Plot of mean volume fraction < (j)(t) > against 
time t for two different shaking amplitudes A. The open cir- 
cles correspond to A = 0.05, while the crosses correspond to 
A = 0.5. 

The mean volume fraction, < (/)(t) >, is evaluated over 
the central 80 per cent of the vertical extent of the particle 
assembly. This central zone is defined by two planes, with 
heights ho = 0.2Z, /ii = 1.8Z, placed symmetrically about 
the (volume weighted) mean particle height Z(t). The gran- 
ular medium is in a box with a hard base, which has periodic 
boundaries in the lateral directions. While all our simulations 
are initiated from stable loose packings with ^ 0.54, we 
note that the rate of structural relaxation increases with the 
shaking amplitude; the time required to reach the steady state 
decreases with increasing A. More importantly however, we 
note that the bulk of the relaxation to the higher packing frac- 
tions occurs for the lower shaking amplitude, which is a direct 
consequence of the predominance of the slow dynamics of 
cooperative relaxation. In this dynamical regime, collective 
structures remain long-standing (their autocorrelation func- 
tions decay only very slowly), and it is possible to show via 
the use of displacement-displacement correlation functions 
that particles move together in dynamical clusters; the en- 
hanced compaction occurs when the net void space trapped 
in these clusters is shaved off. 
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FIG. 2. Plot of mean particle height Z{t) against time t for 
two different shaking amplitudes A. The open circles corre- 
spond to A = 0.05, while the crosses correspond to A = 0.5. 

In Fig. ^, we plot the mean particle height Z(t) against 
time t; this is an important indicator of consolidation as it 
shows the average distance of the particles from the base. 
We notice that the final value of < Z{t) > is smaller for 
smaller shaking amplitudes. In other words, the particles are 
more efficiently settled or 'sedimented' for gentler shaking, 
and consoUdation is thus stronger. 

These two indicators are independent measures of consol- 
idation, but reveal nothing about its nature. Their product 
< Z(t)(/)(t) >, however, provides an insight into whether the 
consolidation is homogeneous or inhomogeneous; this is plot- 
ted in Figure ^. 
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The position Z(t) corresponding to the mean particle height 
at time t is: 
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FIG. 3. Plot of mean (normalised) product of particle 
height < Z(t)(/)(t) > / < Z(oo)(/)(oo) > against time t for 
two different shaking amplitudes A. The open circles corre- 
spond to A = 0.05, while the crosses correspond to A = 0.5. 



Z{t) = [((/)max(t) - (l)o(t))Zmax(t)^ 

+ 0o(t)^o(t)']/20o(O)zo(O) 



(2) 



The position of the consolidation wavefront Zmax{t) rela- 
tive to the central zone, i.e. the region between planes 
/io(t) = 0.2Z(t) and hi(t) = 0.8Z(t), defines the measured 
mean volume fraction < (j){t) > at time t: 

< 0(t) > = 0(0) ^max(t) < ho(t) (3) 

< (/)(t) > — [((/)macc — 0o)2;max(t) 

+ (1.80(0) - O.20^a,)Z(t))]/1.6Z(t), 

9/lo(t) > Zmax{t) > ho{t) (4) 

mace Zmax {t)>9ho{t) (5) 



In the case of homogeneous consoUdation, the volume frac- 
tion and mean particle height are related by a geometrical 
(centroid) condition such that (j){t)Z{t) — constant. However 
our plots of < (j){t)Z{t) > (normalised to the value at t = oo) 
(Fig. ^ show systematic deviations from unity at early times 
and therefore indicate inhomogeneous consolidation at these 
times. 

The details of our findings can be explained by the results 
of a simple theory, which we develop in the next subsection. 



Equations (2) and (3-5) are in qualitative agreement with 
Figures Q and ^ 

In order to explain the non-trivial behaviour of the product 
< (/)(t)Z(t) >, we note that at large times, t ^ oo, we must 
have: 

Zmax{t) = Zo{t) = 2Z(t) (6) 
0maxZ(t) = 0(O)Z(O) (7) 



B. Theory: consolidation waves 

Consider the simplest form of inhomogeneous consolida- 
tion, i.e. the volume fraction changes discontinuously from 
an initial uniform density, 00, to a maximum density, (j)max, 
at each layer of a (three-dimensional) vibrated granular bed, 
following the passage of a 'consolidation wave'. We visualise 
the process of consolidation to be analogous to that which is 
known to occur in sedimentation when particles settle 
onto a sediment at the bottom of a particulate suspension. In 
this way, we imagine that shaken grains settle to the bottom 
of the vibrated bed in a more compacted state than they were 
previously, and that this settling defines the position of the 
front of the consolidation wave. 

Let Zmax (t) be the position of the highest layer with density 
0max at time t; this then defines the position of the front of 
the consolidation wave. Let zo{t) be the position of the upper 
surface (characterised by a density 0o(t) at time t. Assuming 
that consolidation is uniform in the plane perpendicular to the 
direction of vibration (which we choose to be the z-direction) 
, while it is inhomogeneous in the direction of vibration, con- 
servation of mass gives 

(l)max(t)Zmax(t) + 0o(t)(zo(t) - Zmax(t)) 

= 0o(O)^o(O) (1) 



since the consolidation wavefront will have reached the up- 
per surface at large times. Then, using the above equations, 
and simplifying, we find that for short times, the normalized 
product < (j)(t)Z(t) > (which approaches one as t ^ oo), is 
given by 

[(t>{t)Z{t)]/[(l)maxZ{00)] = Z{t)/Z{0) 

= zo{t)/zo{0) + 

{{m/i^max - 0O))( [ (^O(O) - Zo{t))/Zo{0)f) (8) 



The first term on the right- hand-side is expected to be less 
than one in the case of consolidation, since the height of the 
upper surface should decrease with time as the granular bed 
compacts. The second term, in this case, is also clearly pos- 
itive; its size, however, depends on the size of the coefficient 
(j){0) / {(j)max — 0(0)). We have seen in Figure 1 that small val- 
ues of 0max — 0(0) correspond to large shaking amplitudes 
A, and conversely. This results implies, in agreement with 
earlier work on compaction that the overall consolidation 
is larger for lower shaking amplitudes. Thus, for large A the 
normalised product < (j){t)Z{t) > can be greater than one ; for 
small shaking amplitudes, on the other hand, the value of the 
product can be less than one (the value of 0(0) / (0macc — 0(0)) 
is small). These results obtained in the short-time limit of 
this very simple model, are in excellent agreement with the 
results of the simulation (Figure 0) . 
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IV. IRREVERSIBLE AND REVERSIBLE 
BRANCHES IN ANNEALED SHAKING: MORE 
CONSOLIDATION WAVES 

Recent explorations of compaction in a shaken powder have 
involved the analogue of 'annealed cooling' a granular bed, 
prepared at very low density, is subjected to ttap taps at in- 
tensity F; the intensity is then incremented in steps of ST. 
This results in a compaction curve with two branches; an irre- 
versible branch, where increasing the tapping intensity leads 
to increasing packing fractions, and a (so-called) reversible 
branch, where the reverse is the case. We have simulated the 
'annealed shaking' of a granular bed of 1300 particles for 1425 
shakes, such that 6r = 1/500 and ttap = 1, corresponding to 
a continuous variation of shaking intensity. As before, there is 
a hard base in the lower z direction, while there are periodic 
boundary conditions in the x and y directions. The intensity 
is increased from A = 0.05 to A = 1.0, decreased to A = 0.05 
and increased back again to A = 1.0. 

Our results (Fig. show that the packing fraction first 
increases from 0.55 along the irreversible branch before de- 
creasing at high values of A; it then increases as A is reduced 
along the socalled reversible branch. Any further changes of 
shaking intensity cause it to remain on this branch, except for 
a small hysteresis loop; the presence of this hysteresis loop in- 
dicates that even the socalled reversible branch is irreversible 
in reality. We expect that the size of this hysteresis loop would 
increase for large ^F and small ttap, and conversely; the irre- 
versibility of the 'reversible' branch consequently mirrors the 
degree of irreversibility of the annealed shaking process. 



branch, and the last three to A increasing along the reversible 
branch. We note, especially along the 'reversible' trajectories 
that a 'crystalline '-like ordering appears to spread from the 
base as the shake intensity is continuously decreased, which 
has a wavelike character. The presence of this 'consolidation 
wave' mirrors a recent experimental result, where ordering in 
an assembly of rods is seen to originate at the walls and to 
propagates like a wavefront . 
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FIG. 5. Configurations for selected evenly spaced shaking 
amplitudes A, a) on the irreversible branch A = 0.3, 0.6, 0.9, 
b) the 'reversible' branch for A decreasing A = 0.8, 0.5, 0.2, 
and finally c) the 'reversible' branch for A increasing 
A = 0.15,0.45,0.75 



FIG. 4. Plot of packing fraction (f) vs shaking intensity A 
for the annealed shaking of a granular bed. Note the hys- 
teresis loop in the 'reversible branch', showing a fundamental 
irreversibility. 

In order to look more microscopically at what occurs along 
both branches of the compaction curve, we present in Fig. 5 
snapshots of the configurations of the spheres. Each snapshot 
corresponds to a distinct value of the shake intensity A as 
the compaction curves of Fig. ^ are swept. The first three 
configurations correspond to A increasing on the irreversible 
branch, the second three, to A decreasing on the reversible 



Investigations are under way to analyse these structures 
further, with the use of appropriate order parameters |20| in 
the context of an observed crystallisation transition in shaken 
granular beds [EgI. 



V. DISCUSSION 

We have presented a simple model of short time inhomo- 
geneous consolidation, based on the idea of a consolidation 
wave; its analytical solution is in good agreement with our 
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Monte Carlo simulations of vibrated (three-dimensional) par- 
ticle packings. Despite the simplicity of our model, it mani- 
fests an important signature of inhomogeneous consolidation 
for early and intermediate times, in agreement with the sim- 
ulations. Interestingly these are both in agreement with ex- 
perimental data [^, where it was shown that the larger the 
intensity of vibration, the more homogeneous the consolida- 
tion ; notice that the data for the larger vibrational intensity 
asymptote much more quickly (Figure^. Investigations are 
in progress to look in greater detail at the local structures in 
the particle packings, and at the detailed spatial variations 
inherent in the inhomogeneous packings [^. Ongoing work 
also includes a detailed investigation of inhomogeneous den- 
sity fluctuations in a vibrated powder, and their repercussions 
for its power spectrum 

Additionally, we have presented ancillary evidence for con- 
solidation waves via the study of particulate configurations in 
the annealed shaking of a granular bed. In addition to the as- 
pects on which we have reported here, an analytical study of 
a simple lattice-based model of a shaken granular bed is 
expected to yield valuable insights on the precise influence of 
spatial inhomogeneities along the reversible and irreversible 
branches of the compaction curves. 
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